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ABSTRACT 



Context. Red supergiants are massive evolved stars that contribute extensively to the chemical enrichment of our Galaxy. It has 
been shown that convection in those stars gives rise to large granules that cause surface inhomogeneities and shock waves in the 
photosphere. The understanding of their dynamics is crucial to unveil the unknown mass-loss mechanism, their chemical composition 
and stellar parameters. 

Alms. We present a new generation of red supergiants simulations with a more sophisticated opacity treatment done with 3D radiative- 
hydrodynamics C05B0LD. 

Methods. In the code, the coupled equations of compressible hydrodynamics and non-local radiation transport are solved in the 
presence of a spherical potential. The stellar core is replaced by a special spherical inner boundary condition, where the gravitational 
potential is smoothed and the energy production by fusion is mimicked by a simply producing heat corresponding to the stellar 
luminosity. All outer boundaries are transmitting for matter and light. The post-processing radiative transfer code 0PTIM3D is used 
to extract spectroscopic and interferometric observables. 

Results. We show that the relaxation of the assumption of frequency-independent opacities shows a steeper mean thermal gradient in 
the optical thin region that affect strongly the atomic strengths and the spectral energy distribution. Moreover, the weaker temperature 
fluctuations reduce the incertitude on the radius determination with interferometry. We show that ID models of red supergiants must 
include a turbulent velocity calibrated on 3D simulations to obtain the effective surface gravity that mimic the effect of turbulent 
pressure on the stellar atmosphere. We provide an empirical calibration of the ad-hoc micro- and macroturbulence parameters for ID 
models using the 3D simulations: we find that there is not a clear distinction between the different macroturbulent profiles needed in 
ID models to fit 3D synthetic lines. 

Key words, stars: atmospheres - stars: supergiants - hydrodynamics - radiative transfer - Methods: numerical 



1. Introduction 

The dynamical nature of the solar- surface layers, manifested 
for instance in granules and sunspots, has been known for a 
^ long time. With every improvement of ground-based or space- 
c borne instruments the complexity of the observed processes in- 
creased. Red supergiant (RSG) stars are among the largest stars 
in the universe and the brightest in the optical and near in- 
frared. They are massive stars with masses between roughly 
10 and 25 M© with effective temperatures ranging from 3 450 
to 4100K, luminosit ies of 2 0000 to 300 000 Lq, and radii up 
to 15OO7?0 (Levesq ue et al.l [2005). These stars exhibit varia- 
tions in integrated brightness, surface features, and the depths, 
shapes, and Doppler shifts of spectral lines; as a consequence, 
stellar parameters and abundances are difficult to determine. 
Progress has been done using ID hydrostatic models revis- 
ing the Teff- scale an d reddening both at solar and Magellanic 
Clouds metallicities dLevesque et al. 1120051 120061: iMassev et al.] 
120071: iLevesque et al.ll2007l: lLevesqudl2010h but problems still 



remain, e.g. the blue-UV excess that may be due to scatter- 
ing by circumstellar dust or to an insufficiency in the mod- 
els, and the visual-i nfrared effective temperature mismatch 
(iLevesque et al. I [2006b . Finally, RSGs eject massive amounts of 
mass back to the interstellar medium with an unidentified pro- 
cess that may be r elated to acou stic waves and radiation pres- 
sure on molecules (iJosselin & Pl ez 2007), or to the dissipation 
of Alfven waves from magnetic field, re cently discovered on 



RSGs (lAuriereetal.1 [20101: lGrunhutetal.1 120101). as ea r ly sug - 



gested bvlkartmann & AvrettI (11984 : iPiipers & Heard (1 19891) : 
ICuntzl(ll9m 

The dynamical convective pattern of RSGs is then crucial 
for the understanding of the physics of these stars that contribute 
extensively to the chemical enrichment of the Galaxy. There is 
a number of multi wavelength imaging examples of RSGs (e.g. 
a Ori) because of their high lumin osity and large angula r di- 
ameter. Concerning a Ori, Buscher et al.[ ([1990) ; Wilson et al 
1992h: [Tuthill etalJ ([I997h : [Wilson et al.[ ([19971) : [Young et al 



(I2000L 120041) detected the presence of time-variable inhomo- 
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geneities on his surface with WHT and COAST; iHaubois et alJ 
(I2009h reported a reconstructed image in th e H band with 
two large spots; ^hnakaetal. (2009, l201lh attributed mo- 
tions detected by V LTI/ AMBER observations to convection; 
iKervella et alJ (l2009l) resolved a Ori using diffraction-limited 
adaptive optics in the near-infrared and found an asymmet- 
ric envelope around the sta r with a b right plume extending 
in the southwestern region. iTuthill et al . ( 199 7) the presence 
of bright spots on the surface of the s upergia nts a Her and 
a SCO using WHT and lChiavassa et alJ (l2010cl) on VX Sgr us- 
ing VLTI/AMBER. 

The effects of convection and non-radial waves can be rep- 
resented by numerical multi-dimensional time-dependent radi- 
ation hydrodynamics (RHD) simulations with realistic input 
physics. Three-dimensional radiative hydrodynamics simula- 
tions are no longer restricted to th e Sun (for a review on the Sun 
models see iNordlund etaDl2009l) but co ver a substantial por - 
tion of the Hertzsprung-Russell diagram (iLudwig et al.ll2009h . 
Moreover, they have been already extensively employed to study 
the effects of photospheric inhomogeneities and velocity fields 
on the formation of spectral lines and on interferometric observ- 
ables in a n umber of cases, incl i iding the Sun, dwarfs and sub 
giants (e.g. lAsplund et al.lll999l:lAsplund & Garcia Perezl[2001 



Asplund et aP l2009l: ICaffau et al.l l201Ql: iBehara et al.ri201Q : 
Sbordone et al. '2010'), red giants (e.g. 'Chiavas sa et a'lTl2010a : 



Collet et al.l 2007. 2009; Wen de eta l. 2009), asymptotic giant 
branch stars (IChi avassa et al. [2010d) . and red supergiant stars 



(IChiavassa et al ] f2009L 1201051201 Ih . In particular, the presence 
and the characterizati on of the size of convect ive cells on a Ori 
has been showed by IChiavassa et al. 1 (l20T0bh by comparing a 
large set of interferometric observations ranging from the optical 
to the infrared wavelengths. 

This paper is the fourth in a series aimed to present the 
numerical simulations of red supergiant stars with CO^BOLD 
and to introduce the new generation of RSG simulations with 
a more sophisticated opacity treatment. Such simulations are 
of great importance for an accurate quantitative analysis of ob- 
served data. 

2. 3D numerical simulations with CO^BOLD 

2.1. Basic equations 

The code solves the coupled equations of compressible hydro- 
dynamics and non-local radiation transport 



dp 
dt 
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They describe the inviscid flow of density p, momentum pv, 
and total energy eikg (including internal, kinetic, and gravita- 
tional potential energy). Further quantities are the velocity vector 
V, pressure P, and radiative energy flux F . The latter is com- 
puted from the frequency-integrated intensity / for the gray treat- 
ment of opacity and for each wavelength group in the non-gray 
approach (Sect. 12.31) . 

The gravitational potential is spherical. 



where Mpot is the mass of the star to be modeled (see Fig. [3l 
bottom right panel); ro and n are free smoothing parameters: 

; while when r ^ 0, O ^ O = 



when <^ r"^ «c r p O = - 



-^; and for r 



oo, O 



o = - 



. Typically, ro ~ 0.2 



and ri ^ 1.2 for RSG simulations (iFrevtag et al.ll2002h . 



2.2. The code 

The numerical simulations described here are performed with 
CO^BOLD (conservative COde for the Computation of 
compressible COnvectio n in a BOx of L Dimensions, L = 2, 3). 
It uses operator splitting (IStrangI [79681) to separate the various 
(explicit) operators: the hydrodynamics, the optional tensor vis- 
cosity, and the radiation transport. 

The hydrodynamics module is based on a finite volume ap- 
proach and relies on directional splitting to reduce the 2D or 
3D problem to one dimension In the ID steps an approximate 
Riemann solver of Roe-type (lRodll986b is applied, modified 
to account for a realistic equation of state, a non-equidistant 
Cartesian grid, and the presence of source terms due to an ex- 
ternal gravity field. In addition to the stabilizing mechanism in- 
herent in an upwind- scheme with a monotonic reconstruction 
method (typically a piecewise-linear van Leer interpolation), a 
2D or 3D tensor viscosity can be activated. This step eliminates 
certain errors of Godunov-type m ethods dealin g with strong ve- 
locity fields aligned with the grid ('OuirlJ^ 19941) . 

The equation of state uses pre-tabulated values as functions 
of density and internal energy (p, Ci P,Ti, T, s). It accounts for 
HI, HII, H2, Hel, Hell, Helll and a representative metal for any 
prescribed chemical composition. The equation of state does not 
account for the ionization states of metals, but it uses only one 
neutral element to achieve the appropriate atomic weight (in the 
neutral case) for a given composition. Two different geometries 
can be used with CO^BOLD, that are characterized by different 
gravitational potentials, boundary conditions, and modules for 
the radiation transport: 

- The BOX-IN- A- STAR setup is used to model a statistically rep- 
resentative volume of the stellar atmosphere with a constant 
gravitation, where the lateral boundaries are periodic, and 
the radiation transport module relies on a Feautrier scheme 
applied to a system o f long rays (Freytag et all (l2002h. 
Wedemever et al.ll2004 and for analysis, e.g. ICaffau et al.l 



2010h . 



-1/4 



(4) 



- The STAR-IN-A-Box setup is used to model RSG stars of this 
work. The computational domain is a cubic grid equidistant 
in all directions, and the same open boundary condition is 
employed for all side of the computational box. 

As the outer boundaries are usually either hit at some an- 
gle by an outgoing Shockwave, or let material fall back (mostly 
with supersonic velocities), there is not much point in tuning 
the formulation for an optimum transmission of small- amplitude 
waves. Instead, a simple and stable prescription, that lets the 
shocks pass, is sufficient. It is implemented by filling typically 
two layers of ghost cells where the velocity components and the 
internal energy are kept constant. The density is assumed to de- 
crease exponentially in the ghost layers, with a scale height set 
to a controllable fraction of the local hydrostatic pressure scale 
height. The control parameter allows to account for the fact that 
the turbulent pressure plays a significant role for the average 
pressure stratification. The acceleration due to gravity is derived 
from Eq.m Within a radius ro the potential is smoothed (Eq.|4]). 
In this sphere a source term to the internal energy provides the 
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Fig. 1. Gray intensity on one side of the computational cube from the initial sequence of the model St35gm00n05 in Table [TJ The 
axes are in solar radii. The artifacts due to the mismatch between the spherical object and the Cartesian grid become less evident 
with time passing. 



stellar luminosity. Motions in the core are damped by a drag 
force to suppress dipolar oscillations. The hydrodynamics and 
the radiation transport scheme ignore the core completely and 
integrate right through it. 

The code is parallelized with Open Multi-Processing 
(OpenMP) directives. 

2.3. Radiation transport and opacities tables 

The radiation transport step for RSGs' simulations uses a short- 
characteristics method. To account for the short radiative time 
scale several (typically 6 to 9) radiative sub-steps are performed 
per global step. Each sub-step uses only 3 rays (e.g. in the di- 
rections (1,1,0), (1,-1,0), (0,0,1) or (1,0,1), (1,0,-1), (0,1,0) or 
(0,1,1), (0,1,-1), (1,0,0) or 4 rays (along the space diagonals 
(±1,±1,±1)). The different orientation sets are used cyclically. 
The radiation transport is solved for a given direction for both 
orientations. More irregular directions (and more rays) are pos- 
sible but avoided in the models presented here to save computa- 
tional time. The radiation transport operator is constructed to be 
stable in the presence of strong jumps in opacity and temperature 
(even on a coarse grid). 

The frequency dependance of the radiation field in the 
CO^BOLD models can be calculated with two approaches: 

- the gray approximation, that completely ignores the fre- 
quency dependence, is justified only in the stellar inte- 
rior and it is inaccurate in the optically thin layers. The 
Rosseland mean opacities are calculated as a function of 
pressure and temperature (T,P a:ross) and available in 
a 2D table. The necessary values are found by interpola- 
tion in a 2D table. It has bee n merged at around 12 OOOK 
from high-temperature OPAL (Iglesias et al. 1992) data and 
low-temperature PHOENIX fflauschildt et al. 1997i) data by 
Hans-Giinter Ludwig. The models which use theses opaci- 
ties are reported in Table [T] 

- the more elaborate scheme accounting for non-g ray eff'ects 
which is based on the idea of op acity binning (iNordlundl 
119821: iNordlund & DravinsI 119901) . The basi c approxima- 
tion i s the so called mu lti-group scheme (iLudwig et al.l 
119941: IVogler et all l2004 . In this scheme, the frequen- 
cies which reach monochromatic optical depth unity within 
a certain depth range of the model atmosphere will be 
put into one frequency group. The opacity table used in 
this work is sorted in 5 wavelengths groups accordingly 
to the run of the monochro matic optical dep th in a cor- 
responding MARCS (Gustaf sson et al.l |2008[) ID model. 
The corresponding logarithmic Rosseland optical depth are 
-Fcx), 0.0, -1.5, -3.0, -4.5, -cx). In each group there is a 



smooth transition from a Rosseland average in the optically 
thick regime to a Planck average in the optically thin regime, 
except for the group representing the largest opacities, where 
the Rosseland average is used throughout. The implemen- 
tation of non-gray opacities tables has been carried out by 
Hans-Giinter Ludwig and we computed one non-gray model 
(Table©. 



3. Modelling RSG stars 

The most important parameters that determine the type of the 
modeled star are: 

- the input luminosity into the core 

- the stellar mass that enters in Eq.|4]for the gravitational po- 
tential 

- the abundances that are used to create the tables for the 
equation-of- state and the opacities 

In addition, there are a number of parameters that influence the 
outcome of a simulation to some extent, i.e. the detailed formula- 
tion of the boundary conditions, the smoothing parameters of the 
potential, the numerical resolution, detailed settings of the hy- 
drodynamics scheme and the additional tensor viscosity, choice 
of ray directions for the radiati on transport, and the time- step 
control. The model presented in iFreytag et al.l (l2002h relies on 
the same assumptions as the current ones. In the meantime, mod- 
ifications of the code (e.g. restrictions to rays only along axes and 
diagonals avoiding some interpolations in the radiation transport 
module) and faster computers allow models with higher reso- 
lution, frequency-dependent opacity tables, and various stellar 
parameters. These new models show a significant increase in the 
number of small convective cells on the stellar surface. 

The initial model is produced starting from a sphere in hy- 
drostatic equilibrium with a small velocity field inherited from a 
previous model with diff'erent stellar parameters (top left panel 
of Fig. [B. The temperature in the photosphere follows a gray 
T (r) relation. In the optically thick layers it is chosen to be adi- 
abatic. The first frame of Fig. [T] taken just after one snapshot, 
displays the limb-darkened surface without any convective sig- 
nature but with some regular patterns due to the coarse numerical 
grid and the correspondingly poor sampling of the sharp temper- 
ature jump at the bottom of the photosphere. The central spot, 
quite evident at the beginning of the simulation, vanishes com- 
pletely when convection becomes strong. A regular pattern of 
small-scale convection cells develops initially and then, as the 
cells merge, the average structure becomes big and the regular- 
ity is lost. The intensity contrast grows. After several years, the 
state is relaxed and the pattern becomes completely irregular. 
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Table 1. Simulations of red supergiant stars used in this work. 



model 


Simulated 


grid 


grid 


Mpot 


•Menv 


L 






logg 


comment 




relaxed 


[grid points] 


[Re] 


[Mo] 


[Mo] 


[Le] 


[K] 


[Re] 


[cgs] 






time [years] 




















st35gm03n07^ 


5.0 


235^ 


8.6 


12 


3.0 


91 932+1400 


3487+12 


830.0+2.0 


-0.335+0.002 


gray 


st35gm03nl3 


7.0 


235^ 


8.6 


12 


3.0 


89 477+857 


3430+8 


846.0+1.1 


-0.354+0.001 


non-gray 


st36gm00n04 


6.4 


255^ 


3.9 


6 


0.4 


24 211+369 


3663+14 


386.2+0.4 


0.023+0.001 


gray 


st36gm00n05 


3.0 


40P 


2.5 


6 


0.4 


24 233+535 


3710+20 


376.7+0.5 


0.047+0.001 


gray, 



highest 
resolution 



^ Used and deeply analyzed in lChiavassa et al.l (120091 . [20 lObh . Note that a sequence of 1.5 stellar years has been added with respect to the paper. 

Notes. All models have solar chemical composition. The table reports: the simulated time used to average the stellar parameters; the numerical 
resolution in grid points and Rq; the mass used to compute the potential in Eq.|4l total mass actually included in the model (the rest is assumed to sit 
in the unresolved stellar core); the luminosity, effective temperature, radius and surface gravity averaged over both spherical shells and simulated 
time (2"^^ column), the errors are one standard deviation fluctuations with respect to the average over time; comment. 



All memory from the initial symmetry is lost. The influence of 
the cubical box and the Cartesian grid is strongest in the initial 
phase onset of convection. At later stages, there is no alignment 
of structures with the grid nor a tendency of the shape of the star 
to become cubical. 

4. Model structure 

Table [T] reports the simulations analyzed in this work. These 
models are the result of intensive calculations on parallel com- 
puters: the oldest simulation, st35gm03n07, used so far in the 
other papers of this series, took at least a year (this time may 
be difl&cult to estimate due to the queuing system and diflTerent 
machines used) on machines with 8-16 CPUs to provide 9 years 
of simulated stellar time; more recently, we managed to com- 
pute the other simulations of Table [T] on a shorter timescale on 
a machine with 32 CPUs. The model St36gm00n04, which is 
slightly larger in numerical resolution, needed about 2 months 
to compute more than 20 stellar years, while St36gm00n05 (the 
highest numerical resolution so far) about 6 months for about 3 
stellar years. Eventually, St35gm03nl3 needed about 8 months 
to compute 10 years of stellar time despite its lower numerical 
resolution. The CPU time needed for non-gray runs does scale 
almost linearly with the number of wavelength groups. 

4.1. Global quantities 

3D simulations start with an initial model that has an estimated 
radius, a certain envelope mass, a certain potential profile, and a 
prescribed luminosity as described above. During the run, the in- 
ternal structure relaxes to something not to far away from the ini- 
tial guess. We needed several trials to get the right initial model 
for the latest set. The average final stellar parameters are deter- 
mined once the simulation has ended (Fig. [2]). For this purpose, 
we adopted the method reported in Chiavassa et al. ( 2009) (here- 
after Paper I). The quantities are averages over spherical shells 
and over time, and the errors are one-sigma fluctuations with re- 
spect to the average over time. We computed the average temper- 
ature and luminosity over spherical shells, T(r), and L(r). Then, 
we search for the radius R for which 

L(R)/(47tR^) = (tT\ (5) 

where cr is the Stefan-Boltzmann constant. Eventually, the ef- 
fective temperature is Teff = TiR). As already pointed out in 
Paper I, the gray model St35gm03n07 (top row of Fig. [2]) shows 



a drift in the first 2 years (-0.5% per year) before stabilisation 
in the last 2.5 years. This drift in also visible in the radius of 
the non-gray model St35gm03nl3 (second row) even though it 
is weaker (less than 0.3%). The luminosity fluctuations are of 
the order of 1% and the temperature variations < 0.3% for both 
simulations. It must be noted that the last snapshot of the gray 
simulation has been used as the initial snapshot for the non-gray 
model. The two hotter models, st36gm00n04 and st36gm00n05, 
are also somehow related because they both use a gray treatment 
of opacities and st36gm00n05 has been started using the last 
snapshot of St36gm00n04. The temperature (0.5%) and luminos- 
ity (2.2%) fluctuations of the higher resolution model are slightly 
larger than the temperature (0.4%) and luminosity (1.5%) fluc- 
tuations of St36gm00n04. None of the simulation shows a drift 
of the radius in the last years of evolution. 

The bottom row in Fig. [2] show the ratio between the tur- 
bulent pressure and the gas pressure (defined as in Sect. 14.21) . 
This quantity shows that in the outer layers, just above the stel- 
lar radius, the turbulent pressure plays a significant role for the 
average pressure stratification and the radial velocities resulting 
from the vigorous convection are supersonic. 

Figure [3] displays some quantities spatially averaged 
over spherical shells for a snapshot of the gray simulation 
St35gm03n07 in Table [T] Figure |4] shows three-dimensional 
views of some quantities for the same snapshot and simulation. 
Radiation is of primary importance for many aspects of con- 
vection and the envelope structure in a RSG. It does not only 
cool the surface to provide a somewhat unsharp outer boundary 
for the convective heat transport. It also contributes significantly 
to the energy transport in the interior (top left panel of Fig. [3]), 
where convection never carries close to 100% of the luminosity. 
In the optically thick part the stratification is slightly far from ra- 
diative equilibrium and the entropy jump from the photosphere 
to the layers below is fairly large (entropy s in Fig.[3]and Fig.©. 
The He II/He III ionization zone is visible in the entropy struc- 
ture as a small minimum near the surface before the normal steep 
entropy decrease. The opacity peak at around T=13 OOOK just 
below the photosphere causes a very steep temperature jump 
(temperature T and opacity k in Fig. [4]) which is very prominent 
on top of upflow regions (opacity k in Fig. [3]). This causes a den- 
sity inversion (density p in Fig. |?] and Fig|8]) which is a suflftcient 
condition of convective instability resolved, the entropy drop oc- 
curs in a very thin layer, while the smearing due to the averag- 
ing procedure over nonstationary up- and downflows leads to the 
large apparent extent of the mean superadiabatic layer. 
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st35gm03n07 



st35gm03n 1 3 



st36gm00n04 



st36gm00n05 




1000 1200 
Radius [Rq] 



1000 1200 
Radius [Re] 



300 400 500 600 700 

Radius [Re] 



300 400 500 600 700 

Radius [Re] 



Fig. 2. Luminosity, temperature, and radius as a function of time for the simulations of Table [T] from left to right columns 
st35gm03n07, st35gm03nl3, st36gm00n04, and st36gm00n05. The bottom panels are the ratio between turbulent pressure and 
gas pressure for different snapshots. The red vertical lines in all the panels is the approximative position of the radius from Table [T] 



The local radiative relaxation time scale in the photosphere 
is much shorter than a typical hydrodynamical time scale (time 
in Fig. [3]). Numerically, the radiative energy exchange terms 
are rather stiff and prevent large deviations from a radiative 
equilibrium stratification. These terms enforce very short time 
steps (about 600 s per individual radiation transport step) and 
are responsible for the relatively high computational cost of 
this type of simulations. Local fluctuations in opacity, temper- 
ature (source function), and heat capacity pose high demands 
on the stability of the radiation transport module. This is true 
to a lesser degree for the hydrodynamics module due to shocks, 
high Mach numbers, and small scale heights. A side effect of the 
steep and significant temperature jump is the increase in pres- 
sure scale height from small photospheric values to values that 
are a considerable fraction of the radius in layers just below the 
photosphere (pressure scale heigh Hp in Fig. [3]). The non-gray 
model has systematically smaller time scales (Fig. [3l bottom 
row), which causes the increased smoothing efficiency of tem- 
perature fluctuations. 



4.2. Temperature and density structures 

The temperature structure for the 3D simulations is displayed 
in Fig. [5] as a function of the optical depth at A = 5 000 A for 
all the simulations used in this work. We have also computed 
classical I D, LTE, plane-parallel, hydrostatic marcs model at- 
mospheres (iGustafsson et al.ll2008l) with identical stellar param- 
eters, input data, and chemical compositions as the 3D simula- 
tions (Table [2|. It should be noted that the MARCS models do 
not have identical input data and numerical scheme as 3D simu- 
lations but they are a good benchmark for a quantitative compar- 
ison. The optical depth scale in Fig. [5] has been computed with 
the radiative transfer code Optim3D (see Paper I, and Sect. 15.2b 
for all the rays parallel to the grid axes with sl x -y positions 
^ (-^center, ^center) ± (G ^ 0.25) (i.e., a Small square at the center of 
the stellar disk), where ^center and ^center are the coordinate of the 
center of one face of the numerical box, G is the number of grid 
point from Table [H and 0.25 has been chosen to consider only 
the central rays of the box. 

The models using a frequency-independent gray treatment 
show larger temperature fluctuations compared to the non-gray 
case (top right pane l in Fig. O, as was already pointed out in 
iLudwig et al I (1 19941) for local hydrodynamical models, because 
the frequency-dependent radiative transfer causes an intensified 
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st35gm03n07 st35gm03n07 st35gm03n07 




200 400 600 800 1000 1200 1400 200 400 600 800 1000 1200 1400 200 400 600 800 1000 1200 1400 

Radius [Rq] Radius [Rq] Radius [Rq] 




Fig. 3. The two rows show some spatially average quantities (i.e., spherical shells) for a snapshot of St35gm03n07 (Table [T]). Red 
vertical line is the location of the radius. Top row, left panel: luminosities as function of radius in solar radii: continuous curve is the 
total luminosity; the dashed curve is the luminosity due to pressure work; the dotted curve is the luminosity of kinetic energy; and 
the dash-dotted curve is the radiative luminosity. Top row, center: opacity. Top row, right: specific entropy. Central row, left panel: 
pressure scale height. Central row, center: characteristic radiative (continuous line) and hydrodynamical (dotted line) time scales. 
Central row, right: gravitational potential Opot. Bottom row, left panel: random columns opacity and radiative timescale for a given 
snapshot of the gray model St35gm03n07 (black) and non-gray model St35gm03nl3 (red). The five coloured curves correspond to 
the opacity groups. The radiative time scale is estimated for temperature fluctuations with spatial scales of the order of the grid size. 
Bottom row, center: characteristic radiative timescale for the gray (black) and non-gray (red) models. Bottom row, right: temperature 
profiles for the St36g00n04 (255^ grid points, black) and St36g00n05 (401^ grid points, red). 



heat exchange of a fluid element with its environment tending to 
reduce the temperature difl'erences. Moreover, the temperature 
structure in the outer layers of the non-gray simulation tends 
to remain significantly cooler than the gray case (up to -200 
hundreds K, Fig. O, close to the radiative equilibrium. The 
thermal gradient is then steeper in the non-gray model and this 
is crucial for the formation of the spectral lines (see Sect. 15.21) . 
This eff'ect has already b een pointed ou t in the literature and 
a more recent example is IVogled (l2004l) who pointed out that 
the non-gray treatment of opacities is mandatory to compare 
solar magneto-convective simulations with Muram code to the 
observations. 

It is also striking to remark that all the gray simulations largely 
diverge from ID MARCS models, while the case of the thermal 
structure of the non-gray model is more complicated. The outer 
layers show a very good agreement with a cool ID MARCS 
model at 3430K (Fig. [3 top right panel). At Igrsooo ~ 1 (i-^., 



where the continuum flux comes from) the mean 3D temperature 
is warmer than the ID: a hot ID MARCS model at 3700K is 
then necessary to have a better agreement but, however, this ID 
profile diverges strongly for Ig T5000 < 1 • 

While turbulent pressure is naturally included in the RHD 
simulations, it is modeled in ID models assuming a parameteri- 
sation as 

Piurh=PpV^ (6) 

where p is the density, Vt the characteristic velocity, and JS 
is approximatively 1 and depends on whether the motions occur 
more or less isotropically. This pressure is measuring the force 
produced by the kinetic movements of the gas, whether due to 
convective or other turbulent gas motions. In ID models, assum- 
ing spherical symmetry, the equation of hydrostatic equilibrium 
is solved for 
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Fig. 4. Logarithm of temperature (top left panel), density (top right), opacity (bottom left) and entropy (bottom right) of a slice 
through the center from a snapshot of the RSG simulation st35gm03n07 in Table [TJ 



Table 2. ID MARCS models used in this work. 



Tef[ 


M 




Corresponding 


Surface 


[K] 


[Me] 


[cgs] 


3D simulation 


gravity 


3490 


12 


-0.35 


st35gm03n07 


8 


3490 


12 


-0.58 


st35gm03n07 




3430 


12 


-0.35 


st35gm03nl3 


8 


3700 


12 


-0.35 


st35gm03nl3 


8 


3430 


12 


-0.65 


st35gm03nl3 


^eff 


3660 


6 


0.02 


st36gm00n04 


8 


3660 


6 


-0.22 


st36gm00n04 


^eff 


3710 


6 


0.05 


st36gm00n05 


8 


3710 


6 


-0.22 


st36gm00n05 


^eff 



Notes. All ID models are spherical and have been computed with solar 
metallicity. The effective surface gravity (geff) is given at the positions 
of the radii from Fig. [7] and using Eq.[8] The last column displays if the 
surface gravity used to compute ID model is geff or 8- For 8 we took the 
same values as 3D simulations of Table [T] 



VP„ + VPtuA + VP, 



rad 



-P8 



(7) 



where g(r) = ™ (being r the radius, M the mass of the 
star and G the Newton's constant of gravity), Pg = %pT/jdj^o\ 
the gas pressure (% being the gas constant, yUmoi the molecular 
weight, and T th e temp eratu re), and Prad the radiative pressure. 
iGustafsson etaP (l2008h and iGustafsson & PlezI (Il992h showed 
that, assuming P^ad = like in 3D simulations, they could mimic 



the turbulent pressure on the models by using models with those 
effects neglected with an adjusted gravity: 



^eff = 8 



1 



1 +yS 



(8) 



EventuallyjGusta fsson et al.l (l2008l) chose to set Vt = for all 
ID models in their grid, ad vicing those who would have liked a 
different choice to use models with a different mass or gravity, 
according to the recipe given in Eq.[8l 

Having all the necessary thermodynamical quantities in 3D sim- 
ulations averaged over spherical shells, we computed geff using 
Eq.[8]with yUmoi = 1-3 (appropriate value for the atmosphere of 
RSGs). Figure [7] shows the behavior of g/gcs for 3D simulations 
of Table [T] We used the effective surface gravity at the the ra- 
dius position (red vertical lines in the figure) to compute new ID 
models (see Table [2]) . Figure [8] shows the comparison of the 3D 
density structures with ID models. The models with the new ef- 
fective gravity shows better agreements with 3D mean profiles 
than the models with Vt set to zero (i.e., geff = g)- The effects on 
the temperature structures of ID models with different surface 
gravity and equal effective temperature is negligible as displayed 
in all the panels of Fig. [5] 

We conclude that ID simulation of red supergiant stars must 
have turbulent velocity not equal to zero, that accordingly to 
Eq. [8l gives the correct value of surface gravity. The effect of 
the turbulent pressure is a lowering of the gravity that can be 
calibrated on 3D simulations accordingly to Fig. [71 
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st35gm03n07 




st35gm03nl3 



-4.0 



4.0 
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3.5 




3.4 
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-2.0 0.0 

log -^sooo 



Fig. 5. Thermal structures of the simulations in Table [T] as a function of the optical depth at /I = 5000A. Darker areas indicate 
temperature values with higher probability. The solid light blue curve is the average temperature, the red dashed line is the ID 
MARCS model profile with surface gravity g (Table [2l) and the green dotted-dashed line is the ID MARCS model with surface 
gravity geff (see text). In top right panel, the 3D mean thermal profile is compared to a cool MARCS model at 3430K for the outer 
layers and a hot one with 3700K (magenta dotted-triple dashed line) for the continuum forming region (see text). 



5. Non-gray versus gray simulations 

We used the 3D pure-LTE radiative transfer code Optim3D de- 
scribed in Paper I to compute spectra and intensity maps from 
the gray simulation st35gm03n07 and the non-gray simulation 
st35gm03n07 (Table[T]). The code takes into account the Doppler 
shifts due to convective motions. The radiative transfer equa- 
tion is solved monochromatically using extinction coefficients 
pre-tabulated as a function of temperature, density, and wave- 
length. The looki ip tables were comp uted for the solar chemi- 
cal compositions (lAsplund et al.ll2006l) using the same extensive 
atomic and mole cular opacity data as th e latest generation of 
MARCS models (iGustafsson et al.ll2008h . We assumed a zero 
micro-turbulence since the velocity fields inherent in 3D mod- 
els are expected to self-consistently and adequately account for 
non-thermal Doppler broadening of spectral lines. 

5.1. From one to three dimensional models: determination of 
micro- and macroturbulence 

RHD simulations provide a self-consistent ab-initio description 
of the non- thermal velocity field generated by convection, shock 



waves and overshoot that are only approximated in ID mod- 
els by empirical calibrations. Thus, the comparison between ID 
and 3D spectra requires the determination of ad-hoc parameters 
like microturbulence (^mic) and macroturbulence (^mac) veloci- 
ties that must be determined for the ID spectra. 
For th is purpose, we used a similar method as ISteff'en et~an 
(l2009h . Standard ID radiative transfer codes usually apply ^mic 
and ^mac broadening on spectral lines isotropically, indepen- 
dently of depth, and identically for both atomic and molecular 
lines. Hence, to have a good representation of the conditions 
throughout the 3D atmosphere, we selected a set of 12 real Fe I 
lines in the 5000-7000 A range (see Table [3) such that: (i) their 
equivalent width on the 3D spectrum range between few mA to 
150mA; and (ii) with two excitation potentials, low (~ 1.0 eV) 
and high (~ 5.0 eV). Using the ID MARCS models with the 
same stellar parameters (plus the effective surface gravity, geff ) of 
the corresponding 3D simulations (see Table [2l), we fi rst derived 
the abundance from each li ne using Turbospectrum (iPlez et al.l 
varez & Plezll 19981 and further improvements by Plez) 
from the 3D equivalent width. For the non-gray model we used 
the ID MARCS model corresponding to the lower temperature 
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Fig. 6. Mean profiles (top) and temperature diff'erence (bottom) 
of the gray and non-gray simulations from Fig.[5l 



(i.e., 3430K) because its thermal structure is more similar to the 
3D mean temperature profile in the outer layers where spectral 
lines form (see top right panel of Fig. O. Note, that we computed 
the 3D spectra using Optim3D with negligible expected diff'er- 
ences (less than 5%; Paper I) with respect to Turbospectrum. 
The microturbulence velocity was then derived by requiring no 
trend of the abundances against the equivalent widths. The error 
on the microturbulence velocity is estimated from the uncertain- 
ties on the null trend. Once the microturbulence has been fixed, 
the macroturbulence velocity was determined by minimizing the 
diff'erence between the 3D and the ID line profiles, and this for 
three profile types: radial-tangential, gaussian and exponential. 
The error on the macroturbulence velocity is calculated from the 
dispersion of the macroturbulence from line to line. The reduced 

is also determined on the best-fitting ID to 3D line profiles. 
Hence, the automatic nature of the determination ensures an ob- 
jective determination. 

Table |4] presents the derived values of the microturbulence 
and macroturbulence respectively for the 3D gray and non gray 
models. The dispersion of microturbulent velocities is reason- 
ably small, so that the depth-independent microturbulence is 
a fairly good approximation. Note also that the microturbu- 
lence velocity for the gray and non-gray model are rather close. 
Although the macroturbulence profile varies from one to an- 
other, the reduced show that all of them give about the same 
value. Note that the micro and macroturbulence standard devia- 
tions from the average velocities is systematically larger in the 
non-gray model. We found that the non-gray spectral lines gen- 
erally show a more complex profile that aff'ect the line fitting 
with ID models causing larger standard deviations. An exam- 
ple for the non-gray model is shown in Fig.|9l the spectral lines 
look shifted and broadened. The line asymmetries are caused by 
the inhomogeneous velocity field emerging from optical depths 
where lines form. The lower temperature in the outer layers of 
non-gray model (Fig. O could cause smaller pressure and den- 
sity scales heights and consequentially a faster density drops and 



of red supergiant stars with CO^BOLD 
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Fig. 7. Ratio between the g and geff from Eq. [8] for simulations 
St35gm03n07 (solid line) and St35gm03nl3 (dotted line) in top 
panel, and St36gm00n04 (solid line) and St36gm00n05 (dotted 
line) in bottom panel. We used the same snapshots of Fig. [5] and 
[HI The vertical red lines are the approximative positions of the 
radii from Table [T] 



stronger shocks. These shocks would complex the shape of the 
non-gray spectral lines. The macroturbulence parameter used for 
the convolution of ID spectra, reproduce only partly the complex 
profile with larger diff'erences in line wings. 

On a wavelength scale of the size of a few spectral lines, the 
diff'erences between 3D and ID models are caused by the veloc- 
ity field in 3D simulations that aff'ect the spectral lines in terms 
of broadening and shift. ID model do not take into account ve- 
locity fields and ad-hoc parameters for micro- and macroturbu- 
lence are used to reproduce this eff'ect. The microturbulence and 
macroturbulence calibration results are pretty similar for gray 
and non-gray simulations even if the 3D- ID corrections of the 
resulting iron abundances are smaller for the non-gray model. 
The microturbulence values obtai ned in this work a re compara- 
ble, albeit a little bit lower, with ICarr et al" I (l2000h who found 
1.7 < ^mic < 3.5 km/s with a diff'erent method based on ob- 
serv ed CO lines. Our va lues are also in general agreement to 
what lSteff'en et al.l (l2009h found for their RHD Sun and Procyon 
simulations (0.8 < ^mic ^ 1-8). 
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Fig. 8. Density structures of the simulations in Table[TJ Darker areas indicate temperature values with higher probability. The colored 
curves have the same meaning as in Fig.O. 
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Fig. 9. Example of CO first overtone lines computed for the 3D non-gray simulation and the corresponding ID MARCS model of 
Tabled The ^mic turbulent velocity used is 1.5 km/s and the ^mac is 6.4 km/s with a radial-tangential profile (see Table [3]) . 
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Wavelength (/xm) Wavelength (yu,m) 

Fig. 10. Top panel: spectral synthesis of TiO band transition A^O - A (y) (top row) and spectral energy distribution (bottom row) 
for the same snapshots of the gray, St35gm03n07, and non-gray, St35gm03nl3, simulations of Fig. [51 3D spectra are compared to 
the corresponding ID MARCS models (Table [2]) with ^mic and ^mac from Table |4] (radial-tangential profile). 



Table 3. Parameters of the Fe I lines used to extract ID mi- 
croturbulence and macroturbulence and equivalent width, W, 
computed for the gray simulation st35gm03n07 and non-gray 
st35gm03nl3ofTable[I] 



Wavelength x lg(^/) Wgray W non-gray 



A 


eV 




mA 


mA 


5853.148 


1.485 


-5.280 


112.7 


104.0 


5912.690 


1.557 


-6.605 


46.1 


39.8 


6024.058 


4.548 


-0.120 


128.4 


129.8 


6103.186 


4.835 


-0.770 


78.6 


81.1 


6214.182 


1.608 


-7.250 


13.9 


14.6 


6364.366 


4.795 


-1.430 


52.7 


57.8 


6569.216 


4.733 


-0.420 


106.7 


107.1 


6581.210 


1.485 


-4.679 


147.3 


140.7 


6604.583 


4835 


-2.473 


11.4 


14.5 


6712.438 


4.988 


-2.160 


13.6 


18.4 


6717.298 


5.067 


-1.956 


13.8 


17.4 


6844.667 


1.557 


-6.097 


78.9 


71.8 



5.2. Spectral energy distribution 

The cooler outer layers encountered in the convection simula- 
tions are expected to have a significant impact on temperature 
sensitive features. This is the case of, in particular, molecular 
lines (e.g., TiO) and strong low-excitation lines of neu- 



Table 4. ID microturbulence (^mic) and macroturbulence (^mac) 
velocities determined to match the 3D spectra characteristics. 





gray model 




non-gray model 




st35gm03n07 


st35gm03nl3 






(km/s 


) 




1.45+0.09 




1.28 + 0.33 


Log6 (Fe)^D -Log^ (Fe)iD 




0.14 




-0.02 


^mac 


profile 


velocity (km/s) 




velocity (km/s) 


radtan 


6.7 + 0.9 


1.15 


6.4 + 2.7 1.17 


gauss 


10.1 + 1.2 


1.19 


9.7 + 3.5 1.14 


exp 


5.5 + 0.9 


1.21 


5.2 + 2.2 1.22 



Notes. The 3D- ID corrections to Fe abundances, derived for the se- 
lected iron lines with the best matching ^mic, are also indicated. 



tral metals, the line formation regions of which are shifted 
outwards because of t he lower photo spheric temperatures 
(iGustafsson et al.l l2008l: iPlez et al.l Il992h . The importance of 
opacity from molecular electronic transition lines of TiO are 
crucial in the spectrum of M-type stars. Figure [TOl (top row) 
shows the synthesis of a TiO band. The strength of the transition 
depends on the mean thermal gradient in the outer layers 
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('^5000 < 1), where TiO has a large effect, which is larger in the 
non-gray model with respect to the gray one Fig. [6j A shallow 
mean thermal gradient weakens the contrast between strong and 
weak lines and this is visible in the molecular band strength 
which is much stronger for the non-gray model. 
Top right panel of Fig. [TOl shows also that the TiO band strength 
of non-gray model is more similar to the cool ID MARCS 
model at 3430K than the hot one at 3700K. This reflects the fact 
that the 3D mean thermal structure in the outer layers is very 
similar to 1D-3430K model (Fig.O top right panel). 
The approximation of gray radiative transfer is justified only in 
the stellar interior and it is inaccurate in the optically thin layers. 
So far, this approximation has been due to the fact that the RHD 
simulations have been constrained by execution time. However, 
with the advent of more powerful computers, the frequency- 
dependent treatment of radiative transfer is now possible and, 
on a wavelength scale of one or more molecular bands, the use 
of this method is big step forward for a quantitative analysis of 
observations. 

The shape of the spectral energy distribution (SED) reflects 
the mean thermal gradient of the simulations. The absolute flux 
plots in the bottom row of Fig. [TOl display that at lower resolu- 
tion two important conclusions can be retrieved: (i) the spectrum 
based on the gray model is largely difl'erent from the ID spec- 
tra; (ii) the non-gray model shows that the SED, compared to the 
ID MARCS model with the hot temperature (3700K), displays 
almost no distinction in infrared region, and weaker in molec- 
ular bands in the visible region and in the near-ult raviolet re- 
gion. Using the prescriptions by lBessell et al.l (1 19981) for the fil- 
ters BVRJHK, we computed expected colors for the gray and 
non-gray models as well as the corresponding MARCS models 
(Table [5]). The non-gray 3D model shows smaller diff'erences to 
the ID model than the gray one. The stronger radiative efficiency 
in the non-gray case forces the temperature stratification to be 
closer to the radiative equilibrium than in the gray case, where 
convection and wave s have a larger impact on it. 
iJosselin et al.l (|200Q|) probed the usefulness of th e index V-^ as 
a temp erature indicator for Galactic RSGs, but iLevesque et al.l 
(120061) . fitting TiO band depths, showed that V - K SLndV -R 
provides systematic higher effective temperatures in Galactic 
and Magellan clouds RSGs (in average -hl70K and -h30K, re- 
spectively) than the s pectral fitting with ID MARCS models. 
iLevesque et al.l (l2006h concluded that the systematic offset was 
probably due to the limitations of static ID models. Using the 
radius definition of Eq. [S] we find that the eff'ective tempera- 
ture from the SED of 3D non-gray spectrum is reff=3700 K. 
However, in the optical the spectrum looks more like a 3500K 
ID model (judging from the TiO bandheads of Fig. [TOl top right 
panel). Assuming that the 3D non-gray model spectrum is close 
to real stars spectra and using V - K, this leads to a Teff that is 
higher by about 200K than the Tq^ from TiO bands. This goes in 
the right direction given what ILevesque et afl (l2006l) found. 

5.3. Interferometry: visibility fluctuations 

IChiavassa et al ] (l2010bL[20TTh showed that in the visible spec- 
tral region, the gray approximation in our modelling may signif- 
icantly affect the intensity contrast of the granulation. Figure [TT] 
(top row) shows the resulting intensity maps computed in a TiO 
electronic transition. The resulting surface pattern is connected, 
throughout the granulation, to dynamical effects. The emerging 
intensity is related to layers where waves and shocks dominate 



together with the variation in opacity through the atmosphere 
and veiling by TiO molecules. 

The intensity fluctuations are linked with the temperature 
inhomogeneities, which are weaker in the thermal structure of 
the non-gray simulation (see Sect. 14.21 and Fig. [5]). We analyzed 
the impact on the interferometric observables using the method 
described in Paper I. We computed visibility curves for 36 dif- 
ferent rotation angles with a step of 5° from the intensity maps 
in Fig. [TTJ In the plots, we introduced a theoretical spatial fre- 
quency scale expressed in units of inverse solar radii (Rq^). The 
conversion between visibilities expressed in the latter scale and 
in the more usual scale is given by: ['']~^ = [Rq^] • d [pc] -214.9. 
The resulting visibility curves are plotted in bottom row of the 
Fig. [TT] together with the one cr visibility fluctuations, F, with re- 
spect to the average value vis (F = cr/vis). The visibility fluctu- 
ations of the non-gray model are lower than the gray model ones 
for spatial frequencies lower than the first null point (approxi- 
mately 0.0075 ^0^), as a consequence the uncertainty on the ra- 
dius determination (the first null point of the visibility curves 
is sensitive to the stellar radius, iHanbury Brown et^ 1 19741) 
is smaller: -10% in gray model versus ~4% in the non-gray 
one. At higher frequencies, the visibility fluctuations are larger 
in the non-gray model between 0.0075 and 0.002 R~^, then the 
non-gray model fluctuations are systematically lower between 
-0.002 and 0.004 R-\ However, it must be noted that after 0.03 
R~^ (corresponding to 33 7?©, i.e., ~4 pixels), the numerical res- 
olution limit is reached and visibility curves can be affected by 
numerical artifacts (Paper I and Sect.©. 

6. Study on the increase of the numerical resolution 
of simulations 

In this work, we presented two simulations with approximatively 
the same stellar parameters but with an increase of the numeri- 
cal resolution. We analyze the effect of increasing the resolu- 
tion from 255^ grid points of simulation st36gm00n04 to 401^ 
grid points of St36gm00n05 (Table [T]). Due to the large com- 
puter time needed for the high-resolution model, the better way 
to proceed is to compute a gray model with the desired stellar 
parameters first, then the numerical resolution can be increased, 
and eventually the frequency-dependent non-gray treatment can 
be applied. With the actual computer power and the use of the 
OpenMP method of parallelization in CO^BOLD, this process 
may take several months. 

We computed intensity maps at 5000A. This spectral region, 
as for the TiO band, is interested by a granulation pattern with 
high intensity contrast. Figure [121 displays the intensity maps 
permeated by a variety of convection-related structures; while 
the large structures are visible in both images, the higher reso- 
lution model resolves much better the details and smaller scale 
structures appe ar within the large on es (as already shown in a 
previous work, IChiavassa et al.ll2QQ6h . Artifacts due to Moire- 
like patterns are show up in the maps of the higher numerical 
resolution simulation. They are composed of thin ripples up to 
40% brighter than the surrounding points and they are already 
been pointed out in Paper I. They are caused by large opac- 
ities, large variations in opacity, and large changes in optical 
depth between 2 successive columns in the numerical box (see 
K in Fig. [3]). To overcome this problem, we computed visibility 
curves from these images and compared the visibility curves for 
one particular projected baseline from the raw image, and af- 
ter applying a median [3x3] smoothing that effectively erases 
the artifacts. The visibility curves in Fig. [121 show that the sig- 
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Table 5. Photometric colors for RHD simulations of Table [T] and for the corresponding ID models of Table [2j A (3D - ID) shows 
the 3D/1D difference for each index. 



Model 


B-V 


V-R 


V-K 


J-H 


J-K 


H-K 


3D gray, st35gm03n07 


1.680 


0.955 


4.751 


0.859 


1.153 


0.293 


ID, Teff = 3490K, logg = -0.58 


1.881 


1.125 


5.070 


0.931 


1.200 


0.269 


A (3D - ID) 


-0.201 


-0.170 


-0.319 


-0.072 


-0.047 


0.024 


3D non-gray, st35gm03nl3 


1.661 


0.886 


4.356 


0.807 


1.064 


0.258 


ID, Teff = 3700K, logg = -0.35 


1.850 


0.962 


4.202 


0.854 


1.080 


0.226 


A (3D - ID) 


-0.189 


-0.076 


0.154 


-0.047 


-0.016 


0.032 



gray model st35gnn03n07 non-gray nnodel st35gnn03nl3 




Fig. 11. Top row: map of the linear intensity in the TiO band transition A^O - X^A (y) of Fig.[TOl The range is [0; 3.5 x 10^] for the 
gray simulation st35gm03n07 of Table[T]and [0; 4.5 x 10^] ergcm"^ s"^ A"^ for the non-gray simulation st35gm03nl3. Bottom row: 
visibility curves obtained from the maps above. The visibilities are computed for 36 different azimuth angles 5° apart. Bottom right 
panel: visibility fluctuations with respect to the average value for the gray (solid line) and non-gray (dashed line). 



nal from the low-resolution st36g00n04 simulation is larger than 
St36g00n05 for spatial frequencies higher than about 0.07 . 
A lower signal at a particular frequency means better spatial res- 
olution: thus very small structures (corresponding to roughly < 
15-20 Rq) are better resolved in st36g00n05. 
Moreover, the visibility curves are mostly affected by the nu- 
merical artifacts for frequencies higher than -0.05 R^ (corre- 
sponding to 25 Rq, i.e., ~5 pixels) for st36g00n04 and 0.09 R^^ 
(corresponding to 1 1 7?©, i.e., ~4 pixels) for st36g00n05. The in- 
crease of numerical resolution reduces only partly the problem 
in the intensity maps and therefore the cosmetic median filter 
must be still applied. This will not affect the visibilities at lower 
frequencies, which are the only ones observable in practice with 
modern interferometers. 

Apart for the problems on the Moire-like patterns of intensity 
maps, the number of convection related surface structures seems 
to visually increase and possible change size passing from 255^ 
to 401^ grid points (Fig.O- In addition to this, the bottom right 
panel of Fig. [3] shows that going from 255^ to 401^ grid points 



the resolution of the temperature jump does not change so much. 
Thus, the numerical resolution limit has not been reached yet. 
The principal limitation to the computation of larger numerical 
simulation is the computer power. A more complete study con- 
cerning the impact of numerical resolution on the number and 
size of granules is necessary and will be addressed in a forth- 
coming paper where a series of simulations with the same stellar 
parameters and increasing numerical resolution will be analyzed. 

7. Conclusions 

We presented the radiation hydrodynamical simulations of red 
supergiants with the code CO^BOLD together with the new gen- 
eration of models based on a more sophisticated opacity treat- 
ment, and higher effective temperature and surface gravity. 

The actual pattern and the mean temperature structure of 
RSGs simulations are influenced: (i) by differences in the rel- 
ative efficiency of convective and radiative energy transport and 
by the efficiency of radiation to smooth out temperature fluc- 
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Fig. 12. Maps of the square root intensity (to better show the 
structure in the higher resolution simulation) of St36gm00n04 
and St36gm00n05 (Table [B at 5000A together with a visi- 
bility curve extracted for a particular angle (solid curve) and 
the the same visibility's angle after applying a [3x3] median 
smoothing (dashed curve). The intensity range in the maps is [0; 
547.7] erg cm"^ s"^ A'^ 
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tuations; (ii) by the optical depth of the upper boundary of the 
convection zone; (iii) by the extent of convective overshoot. 

The main difference between RSG and Sun-like granulation 
(except for geometrical scales) comes from the altered role of 
radiation: it is much more efficient in transporting energy in a 
RSG than in the Sun. This strongly influences the stratification 
and keeps it closer to the case of radiative equilibrium than its 
inefficient counterpart in the deeper layers of the Sun. IMoreover, 
the strong entropy gradient causes a large buoyancy and render 
the convective motions, that compete with the radiative energy 
transport, more violent in a relatively sharper region over which 
the sub-photospheric entropy jump occurs (compared to a typ- 
ical size of a convective element). Eventually, higher velocities 
are accompanied by larger pressure fluctuations and a stronger 
influence of Shockwaves on the photosphere. 

The most important improvement we brought with this work 
is the relaxation of the assumption of gray opacities through 
multigroup frequency-dependent opacities. The non-gray sim- 
ulation shows: (i) a steeper mean thermal gradient in the outer 
layers that afl'ect strongly the molecular and line strengths that 
are deeper in the non-gray case; (ii) that the general shape of 
the spectrum of 3D non-gray simulation are similar to ID model 
while the 3D gray simulation is largely difl'erent. Hence, we con- 
cluded that the frequency-dependent treatment of radiative trans- 
fer with 5 bins represents a good approximation of the more 
complex profile of ID based on -110000 frequency points. 
IVLoreover, the smaller temperature fluctuations of the non-gray 
model, due to the intensified heat exchange of a fluid element 
with its environment, afl'ect the surface intensity contrast and 
consequentially the interferometric observables reducing the un- 
certainty on the stellar radius determination. 

We also showed that ID models of RSGs must include a tur- 
bulent velocity calibrated on 3D simulations to obtain the efl'ec- 
tive surface gravity that mimic the efl'ect of turbulent pressure on 
the stellar atmosphere. 

We provided an calibration of the ad-hoc micro- and macro- 
turbulence parameters for ID models using the 3D simulations, 
which have a self-consistent ab initio description of the non- 
thermal velocity field generated by convection, shock waves 
and overshoot. We found that the microturbulence velocity for 
the gray and non-gray model are rather close and the depth- 
independent microturbulence assumption in ID models is a 
fairly good approximation, even if the 3D- ID corrections of the 
resulting iron abundances are smaller for the non-gray model. 
We also assessed that there is not a clear distinction between the 
diff'erent macro turbulent profiles needed in ID models to fit 3D 
synthetic lines; nevertheless, we noticed that micro and macro- 
turbulence standard deviations on the average velocities are sys- 
tematically larger in the non-gray model, that shows more com- 
plex line profiles than the gray simulation. This may be caused 
by stronger shocks in the outer layers of non-gray model where 
the pressure and density scales heights are smaller. 

While simulations with more sophisticated opacity treatment 
are essential for the forthcoming analysis, models with a better 
numerical resolution are desirable to study the impact of numer- 
ical resolution on the number and size of granules. We provide 
a first investigation showing that the number of convection re- 
lated surface structures seems to increase and change size pass- 
ing from 255^ to 401^ grid points, thus the numerical resolution 
limit has not been reached yet. Future work will focus on a series 
of simulations with the same stellar parameters and increasing 
numerical resolution. 
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